* The original of this script has been downloaded from 
* https://sites.google.com/site/oscarjorda/home/local-projections
* and adjusted to estimate the response of hours worked to different aggregate shocks
* Cacciatore, Gnocchi, Hauser, February 2022, Daniela Hauser
* Modified by James Cabral, November 2022

clear all
graph drop _all
cap drop all

*add your path here:
cd "C:\Users\jamesCabral\OneDrive - DAZ\CGH_Empirical_Replication"

************Select options
*set number of lags
local p 3
* Choose impulse response horizon
local hmax = 11

********************************
* Choose variables of interest *
********************************

local LHSvarlist mkthrs_nonfarm_pcap nonmkt leisure nonmkt_an leisure_an 
local LHSvarnames `" "Market Hours" "Nonmarket Hours" "Leisure Hours" "Nonmarket Hours, Annual Beta" "Leisure Hours, Anuual Beta" "'

local shockvar vix 
local instrument vix_resid_q  
local shockvar_long "the VIX"

*assign macros to the number of variables being looked at
local LHScount: word count `LHSvarlist'
local shockvarscount: word count `shockvars'

***********************************************
* Assemble quarterly data and identify shocks *
***********************************************

clear

*convert the outputted files from the Matlab program to .dta files
foreach xfilename in YoY_q_annual_beta_nonfarm YoY_q_nonfarm {
	foreach varname in leisure nonmkt {
		import delimited "data/synthetic_series/`varname'_`xfilename'.csv"
		gen daten = 1960.75+(_n*0.25)
		save "data/synthetic_series/`varname'_`xfilename'.dta", replace
		clear		
	}
}

*import the quarterly macro data file
*depending on the version of Stata, there may be problems importing the excel file
*we use the .csv file instead (the data are identical)
import delimited "data/quarterly_macro_data.csv"

drop if missing(daten)

*merge the data files
foreach xfilename in YoY_q_annual_beta_nonfarm YoY_q_nonfarm{
	foreach varname in leisure nonmkt {
		merge 1:1 daten using "data/synthetic_series/`varname'_`xfilename'.dta", nogen
		capture erase "data/synthetic_series/`varname'_`xfilename'.dta"
	}
}

*drop observations that are not in our sample
drop if daten >= 2020

gen time = _n
tsset time

*generate log and per capita variables
gen ln_gdp 					= log(gdp_real)
gen ln_gdp_pcap 			= log(gdp_real/pop_workage)
gen ln_mkthrs_nonfarm		= log(mkthrs_nonfarm_pcap)

*keep sample for which we have vix data (1990 onwards)
reg vix l(1/`p').vix l(1/`p').D1.ln_gdp_pcap l(1/`p').D1.ln_mkthrs_nonfarm l(1/`p').oil l(1/`p').ffr
predict vix_resid_q, residuals

save data/DataQuarterlyLP.dta, replace

********************************************************************************

forvalues i=1/`LHScount'{
local LHS: word `i' of `LHSvarlist'
local LHS_long: word `i' of `LHSvarnames'
	
	clear 
	use data/DataQuarterlyLP

	*merge in raw data for vix, macro uncertainty
	rename daten observation_date
	
	*merge in quarterly identified shocks
	tsset time

	gen LHSvariable = log(`LHS')
	egen SD_`shockvar' = sd(`shockvar')
	gen RHSshock = `shockvar'/SD_`shockvar'
	********************************************************************************
	* Generate LHS variables for the Local Projections (in Levels)

	forvalues h = 0/`hmax' {
		gen LHSvariable_`h' = f`h'.LHSvariable
	}
	
	********************************************************************************
	* Run the Local Projections on Aggregate Hours (in Levels)
	eststo clear
	cap drop b u d Years Zero
	gen Years = _n-1 if _n<=`hmax'+1
	gen Zero =  0    if _n<=`hmax'+1
	gen b=0
	gen u=0
	gen d=0
	
	
	forv h = 0/`hmax' {
				
		ivregress 2sls LHSvariable_`h' l(1/`p').RHSshock l(1/`p').LHSvariable (RHSshock = `instrument'), vce(hac nwest)

		* 90% confidence intervals
		replace b = _b[RHSshock]*100 if _n == `h'+1
		replace u = (_b[RHSshock] + 1.645* _se[RHSshock])*100  if _n == `h'+1
		replace d = (_b[RHSshock] - 1.645* _se[RHSshock])*100  if _n == `h'+1

		eststo
	}

	*create a plot
	*for the market hours plot, we add a y axis label
	if "`LHS'" == "mkthrs_nonfarm_pcap" {
		twoway (rarea u d  Years, fcolor(gs13) lcolor(gs13) lw(none) lpattern(solid)) (line b Years, lcolor(blue) lpattern(solid) lwidth(thick)) (line Zero Years, lcolor(black)), legend(off) title("`LHS_long'", color(black) size(medsmall)) xscale(range(0 11)) xlabel(0(4)8) xtitle("Quarter", size(medsmall)) ytitle("% Change", size(medsmall)) graphregion(color(white)) plotregion(color(white) margin(zero))
	}
	else{
		twoway (rarea u d  Years, fcolor(gs13) lcolor(gs13) lw(none) lpattern(solid)) (line b Years, lcolor(blue) lpattern(solid) lwidth(thick)) (line Zero Years, lcolor(black)), legend(off) title("`LHS_long'", color(black) size(medsmall)) xscale(range(0 11)) xlabel(0(4)8) xtitle("Quarter", size(medsmall)) graphregion(color(white)) plotregion(color(white) margin(zero))
	}
	gr rename `LHS'`shockvar', replace
	
*save the IRFs so that we can make graphs later on
	di "`LHS'"
	gen b_`LHS' = b
	gen u_`LHS' = u
	gen d_`LHS' = d
	keep b_`LHS' u_`LHS' d_`LHS' Years
	keep if !missing(Years)
	save data/graph_`LHS'`shockvar'.dta, replace
		
}

***********************************************
* Graphs to Compare Single Beta, Annual Betas *
***********************************************

clear
drop _all
use data/graph_nonmkt`shockvar'.dta
gen Zero = 0

merge 1:1 Years using data/graph_nonmkt_an`shockvar', keepusing(b_nonmkt_an u_nonmkt_an d_nonmkt_an) nogen
merge 1:1 Years using data/graph_leisure`shockvar', keepusing(b_leisure u_leisure d_leisure) nogen
merge 1:1 Years using data/graph_leisure_an`shockvar', keepusing(b_leisure_an u_leisure_an d_leisure_an) nogen

*overlayed graph for nonmarket hours
twoway (rarea u_nonmkt d_nonmkt  Years, fcolor(gs13) lcolor(gs13) lw(none) lpattern(solid))  (line b_nonmkt Years, lcolor(blue) lpattern(solid) lwidth(thick)) (line b_nonmkt_an Years, lcolor(red) lpattern(shortdash) lwidth(thick)) (line u_nonmkt_an Years, lcolor(red) lpattern(dot) lwidth(thick)) (line d_nonmkt_an Years, lcolor(red) lpattern(dot) lwidth(thick)) (line Zero Years, lcolor(black)), legend(order(2 "Constant {&beta}{sup:j}" 3 "Time-Varying {&beta}{sup:j}{sub:t}") position(6) ring(-1) size(medsmall) rows(2) span) title("Nonmarket Hours", color(black) size(medsmall)) xscale(range(0 11)) xlabel(0(4)8) xtitle("Quarter", size(medsmall)) yscale(range(-0.825 1.52828)) ylabel(-.5(0.5)1.5) graphregion(color(white)) plotregion(color(white) margin(zero)) 
		gr rename nonmkt_`shockvar'_ol, replace

*overlayed graph for leisure hours
twoway (rarea u_leisure d_leisure Years, fcolor(gs13) lcolor(gs13) lw(none) lpattern(solid))  (line b_leisure Years, lcolor(blue) lpattern(solid) lwidth(thick)) (line b_leisure_an Years, lcolor(red) lpattern(shortdash) lwidth(thick)) (line u_leisure_an Years, lcolor(red) lpattern(dot) lwidth(thick)) (line d_leisure_an Years, lcolor(red) lpattern(dot) lwidth(thick)) (line Zero Years, lcolor(black)), legend(order(2 "Constant {&beta}{sup:j}" 3 "Time-Varying {&beta}{sup:j}{sub:t}") position(6) ring(-1) size(medsmall) rows(2) span) title("Leisure Hours", color(black) size(medsmall)) xscale(range(0 11)) xlabel(0(4)8) xtitle("Quarter", size(medsmall)) yscale(range(-0.825 1.52828)) ylabel(-.5(0.5)1.5) graphregion(color(white)) plotregion(color(white) margin(zero))
		gr rename leisure_`shockvar'_ol, replace

di "`shockvar'"
*combine the graphs with the local projection for market hours
gr combine mkthrs_nonfarm_pcap`shockvar' nonmkt_`shockvar'_ol leisure_`shockvar'_ol, rows(1) name(comp_`shockvar', replace) graphregion(color(white)) plotregion(color(white)) subtitle("Shock to `shockvar_long', single beta vs. annual betas",  size(medsmall)) 

gr export output/CGH_figure1.pdf, replace


***********************************************************
* Delete the IRF data files to keep the data folder clean *
***********************************************************

foreach LHS in `LHSvarlist'{
	capture erase data/graph_`LHS'`shockvar'.dta	
}

erase data/DataQuarterlyLP.dta

/* THE END */